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Interim  Report  on 

Numerical  Simulation  of  the  Collapse 
of  an  Underwater  Explosion  Bubble 

ONR  Contract  N00014-87-K-0428 

,  v, 

Purpose:  The  goal  of  this  project  "^s  wdemonstrate^the  usefulness  of  front 
tracking  (J.  Glimm,  et.  al.,  Courant  Institute.  NYUy'as  a  tool  to  study  the  fully 
compressible  axially  symmetric  expansion  and  collapse  of  an  underwater  explosion 
bubble  in  two  spatial  dimensions.  In  the  pilot  project  supported  by  this  contract  a 
one-dimensional  version  of  the  front  tracking  code  is  being  implemented  to  compute 
spherically  symmetric  oscillations.  The  oae  dimensional  version  will  be  used  to  gain 
expertise  in  the  problem,  to  pinpoint  areas  of  difficulty,  and  to  determine  the  feasi¬ 
bility  of  using  the  two  dimensional  front  tracking  method  for  this  problem.^ 

Status  to  date  and  conclusions:  We  have  validatedTa  one  dimensional,  tracked 
random  choice  code  adapted  for  this  pilot  project  on  published  computational  and 
experimental  results.  .These  validation  studies  are  discussed  in  the  Appendix.  tAs  a 
result  of  these  validation  studies,^  we  have  identified  problem  areas  that  exist.  >Iden- 
tified  problem  areas^  include  the  necessity  for  realistic  equations  of  state  for  the 
explosion  products,  implementation  of  effective  boundary  conditions  that  would  take 
into  account  the  progress  of  the  primary  shock  wave  into  the  water  without  the 
necessity  of  keeping  it  within  the  computational  domain,  improved  treatment  of  the 
divergence  at  the  bubble  center  by  analytic  means,  and  correct  initiation  of  the 
explosion  remnants. 

We  have  formulated  plans  for  overcoming  each  of  these  problem  areas.  These 
are  discussed  in  the  next  section  of  this  report.  Finally  we  conclude  that  the  two 
dimensional  computation  of  axially-symmetric  bubble  oscillation  is  easily  possible 
for  deep  water  (i.e.  4.66  kbar  to  426  bar  depths)  explosions.  However  simulation 
of  shallower  depths  (i.e.  74.6  bars)  will  require  some  combination  of  the 
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improvements  specified  below  to  achieve  practical  computational  times. 

Future  improvements:  The  comparisons  with  the  SIN  calculations  of  C.  Mader 
reveal  that  the  use  of  a  polytropic  equation  of  state  for  the  gaseous  explosion  pro¬ 
ducts  is  inadequate.  After  discussions  with  C.  Mader,  we  are  currently  implement¬ 
ing  the  HOM  equation  of  state  for  the  reactive  products.  Two  implementations  are 
planned.  The  first  will  use  the  HOM  equation  of  state  to  generate  data  tables  in  the 
format  of  the  SESAME  material  library.  These  data  tables  can  then  be  used  by  the 
Riemann  solver  developed  by  J.  Scheuerman  of  the  Courant  Institute.  This  will 
allow  random  choice  and  front  tracking  calculations.  The  second  approach  will  be 
the  development  of  an  analytic  Riemann  problem  solver  specific  to  the  HOM  equa¬ 
tion  of  state.  This  will  provide  a  more  robust  and  faster  running  code. 

The  correct  simulation  of  the  boundary  conditions  of  the  oscillating  bubble 
problem  requires  that  the  effects  of  the  primary  shock  wave  traveling  through  the 
water  be  taken  into  account  until  it  is  reduced  to  negligible  strength.  In  all  pub¬ 
lished  compressible  computations  to  date,  this  requirement  has  necessitated  keeping 
the  primary  shock  wave  within  the  computational  domain  until  its  strength  is  negligi¬ 
ble,  at  which  time  it  can  be  terminated  at  a  boundary  that  models  far  field  ambient 
water  conditions.  For  shallower  underwater  explosions,  this  distance  can  be  hun¬ 
dreds  of  bubble  radii.  Most  of  this  distance,  from  the  maximum  bubble  radius  to 
the  far  ranging  primary  shock  is  computationally  wasteful  and  would  best  be 
replaced  by  effective  boundary  conditions.  We  plan  to  remove  the  necessity  of 
tracking  the  primary  shock  over  large  distances  by  using  the  ideas  of  T.  Hagstrom 
and  others  (2]  to  achieve  the  required  effective  boundary  conditions. 

Our  comparisons  with  the  published  work  of  Saito  and  Glass  [4]  have  revealed 
the  inadequacies  of  the  numerical  cutoff  used  at  the  center  of  the  bubble  (radius  = 
0),  which  is  the  common  numerical  implementation  to  avoid  the  divergence  of  the 
Euler  equations  in  spherical  coordinates.  The  treatment  at  the  origin,  for  spheri¬ 
cally  and  cylindrically  symmetric  calculations,  affects  the  dynamics  of  the  calcula¬ 
tion,  especially  if  waves  converging  into  the  origin  produce  strong  reflected  signals, 
as  happens  in  these  bubble  explosion  problems.  We  plan  to  investigate  the  imple¬ 
mentation  of  an  analytic  treatment  of  the  solution  near  the  origin  along  the  line  of 
-Noh’s  work. 
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Finally,  correct  initiation  of  the  explosive  remnants  appea  s  important  in  order 
to  achieve  the  final  level  of  agreement  between  calculations  and  experiments  that  is 
desirable.  This  requires  modeling  the  reactive  phase  of  the  explosion.  We  plan  to 
implement  one  the  models  discussed  by  C.  Mader  [3]. 

■Appendix 

Validation  with  current  published  results:  As  a  first  step  in  validating  the 
tracked  random  choice  code,  we  have  compared  results  with  two  other  published 
one-dimensional,  untracked  random  choice  calculations.  The  first  of  these  is  a  cal¬ 
culation  by  J.  Flores  and  M.  Holt  [1]  of  an  underwater  explosion  at  a  depth  of  one 
foot.  The  calculation  is  fully  compressible,  using  the  inviscid  Euler  equations 
modeling  the  gas  bubble  and  the  water.  The  gas  bubble  is  modeled  with  a  polytro¬ 
pic  gas  equation  of  state,  the  water  with  the  separable  energy,  Tait  equation  of 
state.  The  numerical  method  is  random  choice  combined  with  Sod’s  operator  split¬ 
ting  technique  to  include  the  source  terms  of  spherical  geometry.  The  initial  bubble 
radius  is  1/3  foot  and  the  calculation  is  followed  until  the  spherical  primary  shock 
wave  has  traveled  1  foot,  that  is  until  the  shock  wave  teaches  the  water  surface. 

We  performed  the  same  untracked  calculation,  using  instead  a  stiffened 
polytropic  equation  of  state  for  the  water.  Thermodynamic  profiles  (pressure,  den¬ 
sity,  velocity)  produced  by  the  two  simulations  show  only  small  differences,  attribut¬ 
able  to  the  different  equations  of  state  for  water.  The  comparisons  are  shown  in 
Fig.  1. 

A  second  comparison  was  made  against  the  computation  reported  by  T.  Saito 
and  I.  Glass  [4]  who  performed  an  untracked  random  choice  calculation  of  the 
motion  of  a  helium  bubble  in  air  at  unit  atmosphere  pressure.  Both  the  helium  and 
air  were  treated  as  polytropic  gases.  The  calculation  follows  the  leading  shock  in 
the  air  until  it  has  traveled  five  times  the  initial  bubble  radius,  approximately  the 
time  required  for  one  bubble  oscillation.  Results  again  agree  very  well  when  we 
use  the  same  calculational  scheme  as  Saito  and  Glass.  The  pressure  comparisons 
are  shown  in  Fig.  2.  The  Saito-Glass  experiment  is  however  less  demanding  that  of 
Flores-Holt  as  the  sound  speeds  in  air  and  helium  are  much  more  closely  matched 
than  those  in  gas  and  water. 


-  6  - 


Our  investigations  have  revealed  the  importance  of  the  cutoff  at  the  center  of 
the  bubble  (radius  =  0),  which  is  required  to  prevent  numerical  divergences,  a  point 
glossed  over  in  the  Saito-Glass  report.  It  is  a  known  that  the  treatment  at  the  ori¬ 
gin,  for  spherically  and  cylindrically  symmetric  calculations,  affects  the  dynamics  of 
the  calculation,  especially  if  waves  converging  into  the  origin  produce  a  strong 
reflected  signal.  In  these  bubble  explosion  problems,  such  a  signal  develops  during 
the  collapse  phase.  We  have  shown  that  the  speed  of  the  two  important  waves,  the 
initial  shock  in  the  air  and  the  helium-air  interface  (or,  for  the  explosion  underwa¬ 
ter,  the  shock  in  the  water,  and  the  water-gas  interface)  is  affected  by  these 
reflected  waves.  The  smaller  the  cut-off,  the  slower  the  outgoing  shock  wave  moves 
during  the  compressive  phase.  However,  reducing  the  cut-off  (by  which  we  mean 
actually  calculating  closer  to  the  origin)  also  causes  a  reduction  in  timestep  size 
which  greatly  increases  the  computational  time. 

In  a  final  set  of  validation  studies  we  calculated  the  oscillation  of  a  bubble 
caused  by  detonating  0.55  pounds  of  tetryl  at  three  separate  underwater  depths. 
These  ’'ere  compared  to  calculations  reported  by  C.  Mader  [3]  using  the  SIN  code. 
As  the  SIN  calculations  computed  the  detonation  phase  of  the  explosion,  the  intial 
data  used  by  the  random  choice  code  had  to  be  approximated  from  the  explosion 
data.  The  intial  data  used  are; 

bubble  radius  3.27  cm 

bubble  pressure  0.251  Mbar 

bubble  density  1.7  gm/cm^ 

ratio  of  specific  heats  in  bubble  2.93 

ambient  water  pressure  74.6  bars,  462  bars.  4.66  kbars 

ambient  water  density  l.O  gm/cm^ 

In  comparing  the  SIN  calculations  reported  by  Mader  with  the  tracked  random 
choice  code  results,  two  differences  in  the  approximation  methods  that  are  likely  to 
account  for  differences  in  the  results  should  be  noted.  First,  the  SIN  code  uses  a 
more  realistic  equation  of  state  to  model  the  explosive  products,  whereas  the  ran¬ 
dom  choice  computation  models  the  explosive  products  as  a  gamma  law  gas. 
Second,  the  SIN  code  initiates  a  simulation  by  modeling  the  detonation  phase, 
whereas  the  random  choice  method  is  initiated  at  the  completion  of  detonation  as  a 
high-density,  high-pressure,  underwater  gaseous  sphere  with  velocity  everywhere 
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specified  as  zero.  This  assumption  is  of  course  incorrea;  by  the  time  detonation  is 
complete,  velocities  inside  the  gas  sphere  are  nonzero.  In  addition,  the  random 
choice  calculations  assume  complete  combustion  of  the  explosive,  that  is,  that  the 
initial  bubble  density  is  constant  in  the  bubble  and  equal  to  that  of  undetonated 
tetryl.  A  realistic  detonation  model  would  compute  a  more  accurate  density  profile 
for  the  explosive  products.  These  two  distinctions  between  SIN  and  the  random 
choice  code  used  here  should  account  for  most  of  the  differences  appearing  in  the 
results.  In  the  next  phase  of  the  project,  we  plan  to  implement  the  HOM  equation 
of  state  model  for  the  gaseous  explosive  products  either  in  table-lookup  form  or  in 
analytic  form.  This  upgrade  should  improve  agreement  with  the  SIN  code,  and  with 
experimental  results.  At  a  still  later  stage  of  the  project,  a  more  realistic  initiation 
will  be  introduced,  perhaps  by  modeling  detonation  numerically. 

Figure  3  compares  plots  of  bubble  radius  as  function  of  time  for  the  SIN  com¬ 
putations  and  the  tracked  random  choice  computations  for  tetryl  detonated  at  4.66 
kbars  (~  156.000  feet).  The  random  choice  bubble  has  a  larger  maximum  radius 
and  a  shorter  period  of  oscillation,  but  agreement  with  the  SIN  calculations  at  this 
depth  is  good.  In  particular,  bubble  oscillations  are  damped  and  equilibrium  is 
reached  after  the  same  number  of  periods. 

Figure  4  compares  plots  of  primary  shock  pressures  and  gas-water  interface 
pressures  for  tetryl  detonated  at  4.66  kbars.  The  minimum  interface  pressure  calcu¬ 
lated  by  the  tracked  random  choice  code  is  smaller  than  that  of  the  SIN  calculation, 
a  result  which  is  consistent  with  the  bubble  radius  results.  Pressures  at  the  primary 
shock  agree  very  well. 

Figures  5  and  6  compare  simulations  for  tetryl  at  462  bars  (~  15,500  feet). 
Figure  5  is  a  superposition  of  plots  of  bubble  radius  versus  time.  The  bubble-radius 
profile  ts  calculated  by  the  racked  random  choice  code  has  a  greater  maximum  and 
shorter  period  than  that  of  the  SIN  calculation.  Agreement  with  the  SIN  code  result 
is  within  20%.  Figure  6  compares  plots  of  pressure  at  the  primary  shock  and  pres¬ 
sure  at  the  gas-water  interface.  The  pressure  at  the  bubble-radius  maximum  com¬ 
puted  by  the  tracked  random  choice  code  is  very  nearly  an  order  of  magnitude 
smaller  than  that  computed  by  the  SIN  code,  a  consequence  of  modeling  the  explo¬ 
sion  products  as  a  gamma  law  gas. 
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Figures  7  and  8  are  superpositions  of  plots  of  bubble  radius  versus  time,  pri¬ 
mary  shock  pressure  versus  time,  and  gas-water  interface  pressure  versus  time  for 
tetryl  detonated  at  74.6  bars  (~  2,500  feet).  Figure  7  shows  the  bubble  radius  com¬ 
puted  by  the  tracked  random  choice  code  has  a  smaller  maximum  radius  and  shorter 
period  than  the  SIN  code  calculation.  Figure  8  shows  that  the  interface  pressure  cal¬ 
culated  by  the  tracked  random  choice  code  is  nearly  two  orders  of  magnitude 
smaller  than  that  computed  by  the  SIN  code.  Primary  shock  pressures  agree  until 
the  pressure  in  the  tracked  random  choice  calculation  suddenly  decreases  when  the 
shock  wave  leaves  the  computational  domain.  (The  boundary  pressure  is  held  at 
ambient  value.)  The  smallness  of  the  computational  domain  is  the  reason  why  the 
random  choice  solution  underestimates  maximum  bubble  radius  while  in  previous 
cases  it  over  estimated  maximum  bubble  radius.  This  phenomenon  has  been 
observed  in  earlier  test  computations.  The  primary  shock  wave  reaches  the  com- 
purational  boundary  before  the  bubble  reaches  its  first  minimum.  Further,  the 
difference  between  pressure  behind  the  shock  and  the  ambient  water  pressure  is  too 
large  when  the  shock  wave  reaches  the  boundary. 
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Figure  3.  Superposition  of  plots  of  bubble  radius  (cm)  versus  time  (p-sec)  for  0.55  pounds  of 
tetryl  detonated  at  4.66  kbars  (~  156,000  feet  depth)  as  calculated  by  the  SIN  code  up  to  450 
p.sec  and  by  the  random  choice  code  up  to  *  250  p.sec.  The  bubble  in  the  tracked  random 
choice  calculation  has  a  greater  maximum  and  a  shorter  period.  See  the  text  for  a  discussion 
of  the  source  of  the  differences. 


Figure  4.  Superposition  of  plots  of  pressure  versus  time  at  the  primary  shock  and  at  the  gas- 
water  interface  for  0.55  pounds  of  tetryl  detonated  at  4.66  kbars  (*  156,000  feet  depth)  as 
computed  by  the  tracked  random  choice  code  (lighter  curves)  and  by  the  SIN  code.  The  ran¬ 
dom  choice  calculation  underestimates  the  bubble  pressure  at  the  interface  when  the  bubble  is 
at  its  maximum  c.rtent.  The  pressure  profiles  at  the  primary  shock  agree  more  closely.  Sec 
the  text  for  a  discussion  of  the  source  of  the  differences. 


Figure  5.  Superposition  of  plots  of  bubble  radius  (cm)  versus  rime  (tisec)  for  0.55  pounds  of 
retry  1  detonated  at  462  bars  (=  15,000  feet  depth)  as  calculated  by  the  SIN  code  up  to  1500 
(isec  and  by  the  tracked  random  choice  code  up  to  ~  1000  jisec.  The  bubble  in  the  tracked 
random  choice  calculation  has  a  greater  maximum  and  a  shorter  period.  See  the  text  for  a 
discussion  of  the  source  of  the  differences. 


Figure  6.  Superposition  of  plots  of  pressure  (kbars)  versus  time  (msec)  at  the  primary  shock 
and  at  the  gas-water  interface  for  0.55  pounds  of  tctryl  detonated  at  462  bars  (=  15,000  feet 
depth)  as  computed  by  the  tracked  random  choice  code  (lighter  curves)  and  by  the  SIN  code. 
The  tracked  random  choice  calculation  underestimates  the  bubble  pressure  at  the  interface  by 
nearly  an  order  of  magnitude  when  the  bubble  is  at  its  maximum  extent.  Sec  the  text  for  a 
discussion  of  the  source  of  the  differences. 
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Figure  7.  Superposition  of  plots  of  bubble  radius  (cm)  versus  time  (msec)  for  0.55  pounds  of 
tetryl  detonated  at  74.6  bars  (*  2,500  feet  depth)  as  calculated  by  the  SIN  code  up  to  6  msec 
and  by  the  random  choice  code  up  to  *  4.5  msec.  The  bubble  in  the  tracked  random  choice 
calculation  has  a  smaller  maximum  and  a  shorter  period,  unlike  the  calculations  at  the  previ¬ 
ous  two  depths.  This  results  from  the  primary  shock  wave  leaving  the  computational  domain 
before  the  bubble  reaches  maximum  radius,  consistent  with  earlier  results  that  have  demon¬ 
strated  the  effect  of  using  a  domain  which  is  too  short. 


Figure  8.  Superpositioo  of  plots  of  pressure  (kbers)  versus  time  (msec)  at  the  primary  shock 
and  at  the  gas-water  interface  for  0.55  pounds  of  tetryl  detonated  at  74.6  bars  (*  2,500  feet 
depth)  as  computed  by  the  tracked  random  choice  code  (lighter  curves)  and  by  the  SIN  code. 
The  random  choice  calculation  underestimates  the  bubble  pressure  at  the  interface  by  more 
than  an  order  of  magnitude  when  the  bubble  is  at  its  maximum  extent.  The  pressure  profiles 
at  the  primary  shock  agree  fairly  well  until  the  shock  leaves  the  far  boundary  in  the  random 
choice  calculation.  See  the  text  for  a  discussion  of  the  source  of  the  differences. 


